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Abstract 

Mathematical models of stem cell differentiation are commonly 
based upon the concept of subsequent cell fate decisions, each con- 
trolled by a gene regulatory network. These networks exhibit a multi- 
stable behavior and cause the system to switch between qualitatively 
distinct stable steady states. However, the network structure of such 
a switching module is often uncertain, and there is lack of knowledge 
about the exact reaction kinetics. In this paper, we therefore perform 
an elementary study of small networks consisting of three interact- 
ing transcriptional regulators responsible for cell differentiation: We 
^^ investigate which network structures can reproduce a certain multi- 

^\ stable behavior, and how robustly this behavior is realized by each 

^D network. In order to approach these questions, we use a modeling 

r~ — framework which only uses qualitative information about the network, 

(^ yet allows model discrimination as well as to evaluate the robustness 

^~^ of the desired multistability properties. We reveal structural network 

\ I properties which are necessary and sufficient to realize distinct steady 

^ state patterns required for cell differentiation. Our results also show 

that structural and robustness properties of the networks are related 
to each other. 



1 Introduction 

Stem cells and their potential to give rise to multiple cell types have become 
a focus in systems biology and mathematical modeling throughout the last 



years (Peltier and Schaffer, 2010 MacArthur et al. 2009). Starting from 



*Both authors contributed equally to this work. 



a multi-potent stem cell state, they undergo the process of cell differentia- 
tion which is completed when the cell ends up in a mature cell state. The 
eventual goal in stem cell research is to guide the differentiation of stem 
cells toward a desired cell type safely, and to maintain this state robustly. 
Therefore, a promising yet challenging task is to reveal the processes driving 
cell differentiation by mathematical modeling. 

Cell differentiation is commonly viewed and modeled as a sequence of cell 



fate decisions (Foster et al. , 2009). Thereby, the individual decision steps can 



be represented by modules, each driven by a small gene regulatory network 
of "master" transcriptional regulators (TRs) in which feedback motifs play a 



key role ( Xiong and Terrell , 2003 Huang et al. , 2007 Foster et al. 2009 ) . The 



development of models for these modules however is often severely hampered 
by the lack of detailed knowledge about molecular processes, shifting the 
focus toward qualitative models in order to understand the fundamental 
mechanisms. A variety of publications has contributed to the modeling 
of (stem) cell differentiation, covering a broad range of diverse modeling 



approaches, e.g. Roeder and Glauche (2006); Huang et al. (2007); Narula 



et al. (2010); MacArthur et al. (2008). Usually, specific stem cell systems 



such as hematopoietic or mesenchymal stem cells are considered. Yet, similar 
motifs and system properties can be found among them. For example, all 
these approaches have in common that they represent distinct cell types by 
distinct stable steady states, declaring multistability as a universal system 
property. 

One of several aspects within this field concerns the question of stable 
cell states, and minimalistic models that are capable of reproducing exper- 
imentally observed cell states. Although it has been established to view 
cell differentiation as subsequent decisions driven by multistable switching 



modules (Foster et al. , 2009), the network structure of the single modules 



is often still unclear. We therefore perform an elementary study of small 
networks consisting of three TRs. These networks should be able to exhibit 
one progenitor state and two competing differentiated cell types. We also 
compare two different hypotheses about the steady state levels of the TRs 
in the individual cell types. In order to reveal properties of 3-node networks 
that provide candidate models for the differentiation process, we address the 
following questions: 

• Given three cell type-specific TRs, which networks of interactions be- 
tween the TRs can reproduce the required stable steady states? 

• Among the networks that fulfill the stable steady state requirements, 
are there differences in robustness of the required multistability prop- 
erties with respect to perturbations in the interactions? 

As these are very basic questions, this work builds upon a qualitative mod- 
eling framework which needs only very few assumptions about the reaction 



kinetics. In order to model the influence of a TR on its target, monotonous 
activation and inhibition functions are used, but no exact kinetic parameters 
have to be specified. This framework has been previously used for example 
in 



Chaves et al. (2008); Breindl et al. (2010). 



The outline of this paper is as follows: In Section [2| we describe the 
modeling framework, and define several criteria used to classify the investi- 
gated networks. In Section [3j we present and discuss the results of the model 
analysis. Section |4] gives a short summary and concludes with an outlook 
on future investigations. 

2 Methods 

The modeling framework used for this work is briefly recalled first. We then 
specify the stable steady states to be reproduced by each candidate model, 
and introduce some characteristic measures and a robustness measure. 

2.1 Modeling framework 

The framework is based on ordinary differential equations and makes use 
of only very few modeling assumptions which shall be explained next. As 
we are considering networks of TRs, the influence of such a regulator on 
its target gene has to be described mathematically. It is assumed that 
this influence can be described by monotonous activation and inhibition 
functions. The definitions of these functions are given next. 

Definition 1 An activation (inhibition) function is a function v : [0, oo) ^■ 

[0, N) (fi : [0, oo) -^ (0, N]) with iV G M+ and: 

i) V ([i) is continuously dijjerentiable, 

a) 1^(0) = and z^(x) -^ N as x ^ oo 
(7i(0) = A^ and fi{x) — ;■ as x — )■ ooj, 

Hi) ^{x) (^{x)) is monotonously increasing (decreasing). 

Note, that commonly used kinetics such as Michaelis-Menten or Hill 
kinetics are covered by this definition. Furthermore, it is assumed that each 
of the TRs undergoes a linear degradation. With this, the dynamics of a 
network of n TRs are given by 

Xi = -ki-Xi + fi{x), i = l,...,n. (1) 

In this, X £ M" is the vector of concentrations of the TRs and fi{x) models 
the combined influence of all TRs in a network on the TR Xj. In this study 
we only allow multiplicative combinations of activation and inhibition func- 
tions. As generally the parameters ki and the exact shapes of the activation 



and inhibition functions are not known, (IT]) only captures the interaction 
structure. 

As mentioned in the introduction, our goal is to study all possible net- 
works consisting of three nodes Xj, i £ {1,2,3}. The interaction structure 
of such a network is represented by an interaction matrix 



{^ij} 



ije{l,2,3} 



1 <^=^ Xj activates Xi 
-1 <;=^ Xj inhibits Xi 
<;=^ no interaction. 



(2) 



To give an example, the interaction matrix 



encodes the network 





1 1 




A = 


1 -1 
-10 




±^=-ki- xi + z/i(xi) • 7^2(3:3), 


±2 = - k2 ■ X2 + 1^3{X2) ■ fJ^iixs), 


X3=-k3- 


a;3 + At5(a;i), 





in which the parameters ki and the shapes of the activation and inhibition 
functions are not specified. 



2.2 Specifications of stable steady states 

We are specifically interested in the stable steady states of these networks as 
they determine the cell type. However, in this context we will not consider 
stable steady states as individual points but as forward-invariant sets in 
the state space to account for biological variability. By doing so, only high 
and low levels of a TR are distinguished. To be more precise, we assume 
concentrations x^°", xf^^ and x™^^ with < x[°" < xf^"^ < xf^'^^ to be 
known, and a concentration Xj in the interval X],°™ = [0, x'°™] is considered 



as low. Equivalently, a concentration Xj G X 



high 



high 



, xf"™l is considered 



as high. With this, a stable steady state of ([T]) is considered as a forward- 
invariant hyperrectangular set x = li^. x . . . x X^" with li G {high, low}. For 
ease of notation, the n-tuple of labels (/i, . . . , /„) will be used to refer to this 
hyperrectangular set. 

Let us now turn to our system of three TRs that provide a decision 
module in the cell differentiation process. In this setup, three steady states 
are of special interest: The progenitor state A, a differentiated cell type 
B, and a competing differentiated cell type C. Each of these cell types 
i G {A, B, C} is characterized by a specific 3-tuple of high and low TR- 
levels x(*) = {Lf\l2 ,13 y . One can however find several plausible ways 



to characterize the three different cell types via levels of the three TRs. 
We assume the differentiated cell types x' ^x*^ •* to be characterized by 
a high level of the type-specific TR X2 or X3, respectively. Regarding the 
progenitor state x.^ ', we aim to compare two hypotheses ("specifications" 
of stable steady states, in the following for short SSS-specifications) against 
each other, which represent conceptionally distinct mechanisms: 

(SI) xi corresponds to a progenitor factor, maintaining the progenitor state, 
xi is high in the progenitor state x^ •*, whereas it is low in the differen- 
tiated states x(^\x*-'^\ It has to be downregulated in order to achieve 
cell differentiation. This means the required stable steady states are: 



x 



X 



X 



(A) 
(B) 
(C) 



(high, low, low), 
(low, high, low), 
(low, low, high). 



(3) 



(S2) xi corresponds to a differentiation factor, enabling differentiation, xi 
is low in the progenitor state x.^ ', whereas it is high in the differenti- 
ated states x(^), x'-'-^^ It has to be ttpregulated in order to achieve cell 
differentiation. This corresponds to requiring the stable steady states: 



X 



X 



X 



(A) 
(B) 
iC) 



(low, low, low), 
(high, high, low), 
(high, low, high). 



(4) 



Besides the specified stable steady states, we impose a further require- 
ment on the candidate models: In order to investigate "truly interacting" 
nodes rather than assemblages of isolated nodes or subgraphs, we require 
each 3-node network to be weakly connected. 



2.3 Characteristic measures 

Several characteristic measures are introduced to classify networks with re- 
spect to their interactions, including self-loops, and the types of these in- 
teractions (activating / inhibiting). The number of nonzero entries in the 
interaction matrix A, which equals the zero- norm 1 1 A| |o, is used as a measure 
of the network's connectivity. To measure the preponderance of activating 
versus inhibiting interactions, we use the number of positive entries in A mi- 
nus the number of negative entries in A, which we call the positive-negative 
weight. 

The number of sign-inconsistent (sign- consistent) loops is used to mea- 
sure the degree of inconsistencies (consistencies, respectively) arising from 



network interactions. A loop is either a feedback loop from Xi to itself, or a 
feedforward loop from Xi to Xj. Sign-inconsistencies arise when a feedback 
loop has negative sign, or the two paths of a feedforward loop from Xi to Xj 
have opposite signs. Sign-consistent loops are feedback loops with positive 
sign, as well as feedforward loops with the two paths having the same sign. 



Similar motifs have been outlined e.g. in Ma'ayan et al. (2008). 

To measure how strongly the connectivity is distributed among the three 
nodes, the maximum indegree is computed for each network: 



maxindeg{A) := max < N^ aij > , (5) 

|^ie{i,2,3} , 

which is the largest row-sum, thus the oo-norm ||A||oo. 

2.4 Robustness measure 

Since the goal of this study is to compare networks that can reproduce the 
desired cell types, we computed for each network structure a robustness 
measure TZ that quantifies if and how well this structure can generate a set 
of forward- invariant sets as specified in (Sp). We applied the algorithm 



introduced in Breindl et al. (2010) which is summarized next. 

In order to explain the underlying idea of this method, let us use the 
symbol ip to denote both, activation and inhibition functions. These func- 
tions are indexed ifi^k such that i denotes the regulated TR, i.e., Xi, and k 
enumerates the regulators of Xi . The number of regulators of Xi is denoted 
as Qi. 

Also, let us measure a perturbation of a monotonous function ip as the 
^i-norm of the difference of the original monotonous function if and the 
perturbed monotonous function ip^, i.e.. 



\\cp-ipP\\^= / \^{x)-ipP{x)\dx. (6) 

Jo 

Given a set of desired forward-invariant sets x^, z = 1, . . . ,m, the goal 
of the method is to assign a robustness value 7^ to the interaction structure 
itself, i.e., to the matrix A as in ([2]). This robustness value TZ is defined 
such that it has the following three properties: (i) If a value TZ is assigned 
to a system A, there exists a realization of A with monotonous functions 
(Pi^ki such that the desired sets x^, z = l,...,m, are forward-invariant, 
(ii) If no monotonous function is perturbed by more than 7^, i.e., yi,k : 
\\'fi,k — <^f fclli ^ ^1 it can be guaranteed that forward-invariance of the sets 
x^, z = 1, . . . , ?Ti, is maintained, (iii) For every TZ > TZ, there exist perturbed 
functions (^^^ such that forward-invariance of the sets x^, z = 1, . . . ,m, is 

lost, and for at least one index i, k it holds that W^Pi^k — '^f fclli ~ ^- 




^low ™high /^m^x X 



Figure 1: Illustration of a tube for an activation function. 



The procedure to compute TZ for a given network structure A and spec- 
ification x^, z = 1, . . . ,m, is outlined next. To this end, two definitions are 
given next. 

Definition 2 The 3-tuple of pairs of positive real numbers 



((a;^°™. 



-ylow^ 



^^high^^high)^^^. 



max max 

) 7 



)) such that 7^°™ < -/^'^^ < 7" 



and a;'°™ < x^^^^ < x™^"^ is called tube for an activation function. Equiva- 

lently, the 3-tuple of pairs of positive real numbers 

T,, 



((xl°*, 7'^'g'^), {x^'^^, 7^°"), {x"" 



S7™'')) such that 7"^^^^ < 7'°™ < j 



and X < X ^^ < x™^^ is called tube for an inhibition function. 

Definition 3 An activation function v (inhibition function fi) is said to 
satisfy a tube T^, (T^), denoted asv'^T^y (fi \= T^) if the following inequalities 
hold. 



low 



Mx < X 
Vx > x^'s^ 
Vx < x'"^^ 



v{x) < 7'°" (/i(x) > -f^'>i^) 
i/(x) > 7'^'s^ (/i(x) < 7!°™) 
z^(x) < 7°^^=^ (nix) > 7™"^) 



(7) 
(8) 
(9) 



This means that a tube for a monotonous function restricts the space 
where the graph of this function can evolve (see also Figure [I|. Note that the 
x-values of the tubes are given by the specification of the forward-invariant 
sets. It is now possible to formulate conditions on the 7-values of the tubes 
such that the specified sets are forward-invariant if all functions lie inside 
their tubes, i.e. \/i, k : ipi^k 1= T*' , where T can mean Ty or T^ as required, 
and the tubes are indexed as the according monotonous functions. 



Proposition 1 (Breindl et al. (2010)) Given a set :x. = I^J- x 



with li G {low, high} and given tubes T^ that satisfy the conditions 
Vi € {1, . . . n} : —ki • Xj + 7^. 



UA 



7. > 

—T-Si ~ 



X li: 



(10) 



where x,- = min _^i- Xj, and, with Xj denoting the argument of (pi^^, 



XikziJ-x 



Li,k 



if (Pi,k = ^i,k A G Ix. 

min{7 : (x,7) G T*''^ A x G TjA otherwise 



and 

Vi G {1, . . . n} : -h ■ Xi + 7,^1 • . . . • 7^^^^ < (11) 

where Xi = max i- Xi, and, with Xj denoting the argument of ipi^kj 

7i^fc = max{7 : (x, 7) G T^^ A Xj G X^.}. 

/f Vi, A; : Lfi^k 1= T*' , t/ien the set x is forward-invariant for system (Mj. 

Starting from this result, the robustness measure can now be introduced. 
As forward-invariance for the specified sets can be guaranteed as long as 
\/i,k : ifi^k 1= T*' , the goal is to find, for each tube T*' , the function (pi^k^ 
which is best centered to the tube T*''^. In this context, best centered means 
that (pi^k has to be perturbed more than any other function (pi^k ^ T"^'^ 
in order to violate at least one of the constraints ([7]|9]) (denoted as ipi^k ^ 
T*' ). The minimal perturbation of this best centered function (pi^k in order 
to achieve a violation of the tube constraints is given as the result of the 
optimization problem 

^max^^i.fc)^ sup 7^'^'"((^i,fc,r'*=), (12) 

with 

n^^{^.,u,r^^) = inf ||9.,fc-v^fj|i. (13) 

Then, the robustness measure is obtained by maximizing the smallest value 
7^™'^^(T*' ) of all tubes in the network over all tubes T*' that satisfy Propo- 
sition [TJ 

Definition 4 Given a system Tw with unspecified activation and inhibition 
functions. The robustness measure IZ for the system is defined as 

7^=maxmin7^'"^^(r''=) 

T*.* i,k ^YA) 



s.t: Vx"^ : (10) and (11) hold. 



This measure has the interpretation given at the beginning of this section. 
Furthermore, if there exists no realization of A such that the desired sets 



are forward-invariant, the optimization problem (14) is infeasible. It could 



furthermore be shown that for networks where the functions fi{x) are prod- 



ucts of activation and inhibition functions, problem (14) results in a convex 



optimization problem (Breindl et al. , 2010) 
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Figure 2: (Upper row) Distributions of characteristic measures of all 19008 
possible weakly-connected 3- node networks (blue), of the 206 models re- 
producing the SSS-specification (SI) (red, x20), and of the 242 models re- 
producing (S2) (green, x20). (Lower row) Mean robustness values for each 
characteristic measure, over the 206 models reproducing (SI) (red), and over 
the 242 models reproducing (S2) (green). 

3 Results 

3.1 General observations 

Constructing all possible 3-node networks with interactions as specified in 
Q, i.e. each interaction having exactly one value S {—1,0, 1}, yields 3^'^ = 
19683 possible networks. Requiring the networks to be weakly connected 
leaves 19008 networks. Out of these candidate networks, the algorithm 



from Section 2.4 identifies 206 models that meet the specification (SI), i.e. 
1.08% out of all networks; and 242 models that are able to reproduce the 
specification (S2), i.e. 1.27% out of all networks. 

The distributions of the characteristic measures as outlined in Section [2] 
are shown in Fig. ^ Models for (SI) (Fig. ^ upper row, red) and models for 
(S2) (Fig. [2] upper row, green) show similar distributions of nonzero entries. 
The distribution of the positive-negative weight for models of (SI) is shifted 
to the left, i.e. toward more negative entries, whereas for models of (S2) 
it is slightly shifted to the right, i.e. toward more positive entries. These 
differences are in accordance with the number of inconsistent loops which 
require negative entries: Models for (SI) show a distribution of inconsistent 
loops up to high numbers and only intermediate amount of consistent loops, 
whereas models for (S2) tend to low numbers of inconsistencies and more 
consistent loops. Also, models for (SI) (red) are found for all numbers of 
inconsistent loops, whereas models for (S2) (green) are restricted to up to 



8 inconsistent loops. In contrast, the (Sl)-niodels are more restricted to 
certain numbers of consistent loops than the (S2)-models are. 

Summarizing, there are less networks supporting the hypothesis of a pro- 
genitor factor (SI), and they show a tendency toward negative (inhibiting) 
interactions and hence inconsistent loops. The number of networks provid- 
ing models for the hypothesis of a differentiation factor (S2) is higher. On 
average, these networks have about equal number of positive and negative 
entries, thus reducing the amount of inconsistent loops in favor of consistent 
loops. 

3.2 Structural properties 

Regarding the network structure, we found that all models that can repro- 
duce a SSS-specification can be classified leading to sufficient and necessary 
conditions on the interaction structure. These findings are summarized in 
the following observations. 

Observation 1 Let A be an interaction matrix M). // the network repre- 
sented by A is such that each node Xi, i G {1,2,3} together with its incoming 
links is contained in Fig. [^ then there exists a realization of A that can re- 
produce the SSS-specification (SI). Also the other direction holds: For each 
weakly connected 3-node network that can be assembled from nodes xi to- 
gether with its incoming links contained in Fig. [^ there exists a realization 
that can reproduce the SSS-specification (SI). 

An equivalent observation can be made for the network that can repro- 
duce the SSS-specification (S2). 

Observation 2 Let A be an interaction matrix M). // the network repre- 
sented by A is such that the node xi and the nodes Xj, j G {2,3} together 
with their incoming links are contained in Fig. [^ then there exists a real- 
ization of A that can reproduce the SSS-specification (SI). Also the other 
direction holds: For each weakly connected 3-node network that can be as- 
sembled from nodes xi and Xj, j S {2,3}, together with their incoming links 
contained in Fig. [^ there exists a realization that can reproduce the SSS- 
specification (S2). 

Prom Fig. [3] it can also be explained why the positive- negative weight of 
models for SSS-specification (SI) tends toward a larger number of negative 
entries: All interactions between two different nodes in the network need to 
be negative (inhibiting). Only self-loops are allowed to be activating. For 
SSS-specification (S2) on the other side also activating links from xi to the 
other nodes are possible. 

Using these "building blocks" as listed in Fig. [3] it is possible to construct 
all possible 3-node models that can, for an appropriate choice of activation 

10 



(Sl)-models, every node Xi, i G {1,2,3} 

o n o o o 

U^k U^k U^^ ^^B vjf^ Lfiv 

r \ ^ 7^ r\ 

(S2)-models, node xi 

n n n n 
^ ^ *^ ^ 

(S2)-models, nodes Xj, j G {2,3} 



u^ ^^1 — ^^ ^^ 



:,>!_ ^3>|— ^J^[- ^Mh 



Figure 3: Complete enumeration of all "building blocks" . Networks assem- 
bled from these building blocks are able to reproduce the respective SSS- 
specification. (Upper row:) Links entering Xj from bottom left and right 
represent links from nodes Xj, j ^ i and x^, i ^ k ^ j. (Middle row:) Links 
entering xi from below left and right represent links from nodes X2 and x^- 
(Lower row:) Links entering from top right represent links from xi. Links 
entering from right represent links from Xk, k £ {2,3}, k ^ j. 

and inhibition functions, reproduce the SSS-specification (SI) or (S2). To 
see possible applications of our approach, note that a switching module for 
mesenchymal stem cell differentiation based on literature data and analyzed 
in detail in Schittler et al. (2010!) satisfies specification (SI) and can indeed 



be composed from the building blocks presented in Fig. [3| By applying the 
algorithm from Section [2^ it is in principle possible to compute the building 
blocks for any SSS-specification. The algorithm is also directly applicable 
to larger networks. 

3.3 Robustness properties 

The mean robustness value TZ for all characteristic measures is shown in 
Fig. M (lower row). Two important observations can be made: There is a 
clear negative correlation of 7^ with increasing nonzero entries, and a clear 
positive correlation of TZ with increasing positive-negative weight. Both 
observations are true independent of the specific SSS-specifications (red, 
green). In contrast, no clear correlation of TZ versus the number of inconsis- 
tent or consistent loops can be detected. 



11 



It is stated by Kwon et al. (2008) that perturbations on nodes involved in 



fewer (more) than the average number of feedback loops have a lower (higher, 
respectively) impact on the overall network robustness. Thus, we investi- 
gated whether there is a similar relationship between the maxindeg and the 
robustness value TZ. The results of our analysis are shown in Fig. |4j Among 
all (Sl)-reproducing models, only four different robustness values have been 
assigned while for the (S2)-reproducing models seven different robustness 
values have been computed. Furthermore, the (Sl)-models are more abun- 
dant in higher robustness values, whereas the majority of (S2)-models has 
a lower robustness value (cf. Fig. [4]). That is, models for hypothesis (S2) 
require a finer "tuning" than models for (SI), biologically pointing toward 
more vulnerable regulatory networks. 

Yet, despite these differences, the similarities between the two distribu- 
tions are eminent which suggests a strong relation between the network's 
structural properties and their robustness: Independent of the specification 
(SI) or (S2), networks with a higher maxindeg tend to have lower TZ val- 
ues. This means that, the more incoming interactions a node has, the more 
vulnerable the network gets to perturbations of these interaction. This is in 



0.048 

0.042 
0.038 

0.023 



Models for (SI) 


0.048 


Models for (S2) 


n=17« 


n=3 • 


n=37« 


0.042 


n=7 • 


n=61« 


0.038 


n=9 • - 


n=91« 


0.027 
0.025 
0.023 


n=51« n=24« - 
n=39« - 
n=45« - 




0.02 


n=64« - 


1 2 3 

maxindeg 


1 2 3 

maxindeg 



Figure 4: The robustness value TZ versus the maxindeg for all models. 
For each (maxindeg, TZ), there are n models which have these values of 
maxindeg, TZ. (Left) For specification (SI), there are four different robust- 
ness values, each of them uniquely assigned to a maxindeg value. (Right) 
For specification (S2), there are seven different robustness values. The as- 
signment of a given 7^ to a maxindeg value is not unique. 
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accordance with the findings of Kwon et al. (2008) about network structure 
and their notion of robustness. 

4 Conclusion 

The main focus of this work was to analyze which networks consisting of 
three transcriptional regulators are in principle able to reproduce a set of 
required stable steady states, to search for interaction patterns that are 
necessary in order to meet these requirements, and to compute and compare 
the robustness properties of these networks. Still, the investigations raise 
new questions and possible future work, which we point out briefly in this 
concluding outlook. 

Sign-inconsistencies make a system non-monotone, and the more incon- 
sistent loops are contained in a network, the "farther" the system is from 



monotonicity. Ma'ayan et al. (2008) argue that intracellular systems will 



be "close-to-monotone". Kwon et al. (2008) report that a higher number 



of positive feedback loops and a smaller number of negative feedback loops 
result in a higher robustness. These arguments are in accordance with the 
result here that positive (negative) entries increase (decrease) the robust- 
ness of a network. But in contrast to the number of positive / negative 
entries, the number of inconsistent and consistent loops does not show a 
clear effect on the robustness as defined here. What is more, inconsistent 
loops are even required for some SSS-specifications to be met. These studies 
use different notions of robustness, and their interrelation should be worth 
further investigations. 

Furthermore, one might argue that especially in cell differentiation, stim- 
ulus inputs and transitions between stable steady states rather than just the 
existence of stable steady states play an important role. However, although 
these aspects were not considered in this study, the approach in this pa- 
per provides a first selection step of networks with respect to one necessary 
condition (out of several), namely by checking whether the required stable 
steady states can be reproduced. These aspects provide plenty of oppor- 
tunities for separate investigations, and we hope to address these issues in 
future work. 
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